biosignalsnotebooks | download button biosignalsnotebooks | binder server biosignalsnotebooks | home button biosignalsnotebooks | contacts button biosignalsnotebooks | github button

Event Detection - Muscular Activations (EMG)
Difficulty Level:
Tags detect☁emg☁tkeo

Skeletal muscle activation is, in normal conditions, a voluntary process triggered by a nervous impulse that propagates along motor neurons until the desired muscle.

When the nervous impulse reaches sarcolemma (muscle fiber membrane) the depolarisation/repolarisation continues and the changes in membrane potential can be monitored with specialised sensors placed at skin surface.

For contracting a muscle, a large set of motor units needs to be activated, so that the acquired EMG signal is the sum of their elementary potential changes. Because of this "summation" process, EMG seems to be a little "anarchic", and the essence of EMG signal processing is in study the activation zones (bursts).

So, burst detection is an important processing step, which can be achieved by single or double threshold algorithm, generally preceded by a smoothing phase.

In this Jupyter Notebook it will be presented a single threshold algorithm, that includes the Teager-Kaiser Energy Operator (TKEO) in his implementation.

1 - Importation of the needed packages and definition of auxiliary functions

In [1]: 
# biosignalsnotebooks python package
import biosignalsnotebooks as bsnb

# Numpy package is dedicated to simplify the work (operations between) with arrays/lists
from numpy import cumsum, concatenate, zeros, linspace, average, power, absolute, mean, std, max, array, diff, where

# Scientific packages
from scipy.signal import butter, lfilter
from scipy.stats import linregress
In [2]: 
# Base packages used in OpenSignals Tools Notebooks for ploting data
from bokeh.plotting import output_file, show
from bokeh.io import output_notebook
from bokeh.layouts import gridplot
output_notebook(hide_banner=True)

2 - Load of acquired EMG data

In [3]: 
# Load of data.
data, header = bsnb.load_signal("emg_bursts", get_header=True)

3 - Identification of mac address of the device and the channel used during acquisition

In [4]: 
channel = list(data.keys())[0]
In [5]: 
print ("\033[1mChannel: \033[0m " + str(channel))

Channel:  CH3

4 - Storage of sampling rate and acquired data inside variables

In [6]: 
# Sampling rate and acquired data
sr = header["sampling rate"]
device = header["device"]

# Signal Samples
signal = data[channel]
time = bsnb.generate_time(signal)

5 - Binarisation of EMG signal
5.1 - Preprocessing Steps

In [7]: 
# [Baseline Removal]
pre_pro_signal = signal - average(signal)

# [Signal Filtering]
low_cutoff = 10 # Hz
high_cutoff = 300 # Hz

# Application of the signal to the filter.
pre_pro_signal = bsnb.aux_functions._butter_bandpass_filter(pre_pro_signal, low_cutoff, high_cutoff, sr)

5.2 - Application of TKEO operator

\begin{equation} TKEO[i] = \begin{cases} EMG_{original}[i], & \mbox{if } i=0 \mbox{ or } i=N-1 \\ EMG_{original}[i]^2 - (EMG_{original}[i + 1] \times EMG_{original}[i - 1]), & \mbox{otherwise}\end{cases} \end{equation} ... being $N$ the number of acquired samples.

In [8]: 
# [Application of TKEO Operator]
tkeo = []
for i in range(0, len(pre_pro_signal)):
    if i == 0 or i == len(pre_pro_signal) - 1:
        tkeo.append(pre_pro_signal[i])
    else:
        tkeo.append(power(pre_pro_signal[i], 2) - (pre_pro_signal[i + 1] * pre_pro_signal[i - 1]))
In [9]: 
bsnb.plot([list(time), list(time)], [list(signal), list(tkeo)], legend=["Original EMG", "TKEO Signal"], grid_plot=True, grid_lines=1, grid_columns=2, opensignals_style=True, x_axis_label="Time (s)", y_axis_label=["Raw Data", "Raw Data"])

5.3 - Smoothing Phase
5.3.1 - Definition of Constants

In [10]: 
# Smoothing level [Size of sliding window used during the moving average process (a function of sampling frequency)]
smoothing_level_perc = 20 # Percentage.
smoothing_level = int((smoothing_level_perc / 100) * sr)

5.3.2 - Signal Rectification

In [11]: 
# [Signal Rectification]
rect_signal = absolute(tkeo)

5.3.3 - Application of the rectified signal to a first smoothing stage

In [12]: 
# [First Moving Average Filter]
rect_signal = bsnb.aux_functions._moving_average(rect_signal, sr / 10)

5.3.4 - Application of the rectified signal to a second smoothing stage

In [13]: 
# [Second Smoothing Phase]
smooth_signal = []
for i in range(0, len(rect_signal)):
    if smoothing_level < i < len(rect_signal) - smoothing_level:
        smooth_signal.append(mean(rect_signal[i - smoothing_level:i + smoothing_level]))
    else:
        smooth_signal.append(0)
In [14]: 
bsnb.plot([list(time), list(time)], [list(tkeo), list(smooth_signal)], legend=["TKEO Signal", "Smoothed Signal"], grid_plot=True, grid_lines=1, grid_columns=2, opensignals_style=True, x_axis_label="Time (s)", y_axis_label=["Raw Data", "Raw Data"])

5.4 - Definition of the detection threshold

In [15]: 
# [Threshold]
avg_pre_pro_signal = average(pre_pro_signal)
std_pre_pro_signal = std(pre_pro_signal)

Accordingly to the method proposed by Solnik , threshold value can be defined as:

\begin{equation} threshold = \mu_{\scriptsize EMG} + h\sigma_{\scriptsize EMG} \end{equation}

being $\mu_{\scriptsize EMG}$ the average EMG value, $\sigma_{\scriptsize EMG}$ his standard deviation and $h$ a variable that defines the threshold level.

To ensure that threshold level 100% is not bigger than the smooth_signal and level 0 % is not smaller than the smooth_signal we need to define a normalisation regression function.

In [16]: 
# Regression function.
def normReg(thresholdLevel):
    threshold_0_perc_level = (- avg_pre_pro_signal) / float(std_pre_pro_signal)
    threshold_100_perc_level = (max(smooth_signal) - avg_pre_pro_signal) / float(std_pre_pro_signal)
    m, b = linregress([0, 100], [threshold_0_perc_level, threshold_100_perc_level])[:2]
    return m * thresholdLevel + b

Calculation of two threshold values

In [17]: 
# Chosen Threshold Level (Example with two extreme values)
threshold_level = 10 # % Relative to the average value of the smoothed signal
threshold_level_norm_10 = normReg(threshold_level)

threshold_level = 80 # % Relative to the average value of the smoothed signal
threshold_level_norm_80 = normReg(threshold_level)

threshold_10 = avg_pre_pro_signal + threshold_level_norm_10 * std_pre_pro_signal
threshold_80 = avg_pre_pro_signal + threshold_level_norm_80 * std_pre_pro_signal
In [18]: 
fig_list = bsnb.plot([list(time), list(time), list(time)], [list(int((0.5 * max(tkeo)) / max(smooth_signal)) * array(smooth_signal)), list(smooth_signal), list(smooth_signal)], title=["Original and Smoothed Signals", "Threshold 10%", "Threshold 80%"], grid_plot=True, grid_lines=2, grid_columns=2, hor_lines=[[], [threshold_10], [threshold_80]], opensignals_style=True, show_plot=False, x_axis_label="Time (s)", y_axis_label=["Raw Data", "Raw Data", "Raw Data"], get_fig_list=True)
In [19]: 
fig_list[0].line(time, tkeo)
grid_plot = gridplot([[fig_list[0]], [fig_list[1], fig_list[2]]], **bsnb.opensignals_kwargs("gridplot"))

show(grid_plot)

The threshold level of 10 % is chosen for our application, because, as can be seen in the previous figure, none activation period is completely below the threshold line.

5.5 - Binarisation of the smoothed signal

In [20]: 
# Generation of a square wave reflecting the activation and inactivation periods.
binary_signal = []
for i in range(0, len(time)):
    if smooth_signal[i] >= threshold_10:
        binary_signal.append(1)
    else:
        binary_signal.append(0)
In [21]: 
bsnb.plot([list(time), list(time)], [list(tkeo / max(tkeo)), list(binary_signal)], legend=["TKEO Signal", "Binarised Signal"], grid_plot=False, opensignals_style=True, x_axis_label="Time (s)", y_aAxisLabel="Raw Data")

5.6 - Begin and end of activation periods

All upward transitions (0 to 1) define the beginning of an activation period and all downward transitions (1 to 0) establishes his end

In [22]: 
diff_signal = diff(binary_signal)
act_begin = where(diff_signal == 1)[0]
act_end = where(diff_signal == -1)[0]
In [23]: 
print ("\033[1mBurst Begin times: \033[0m\n" + str(array(time)[act_begin]) + " s")
print ("\033[1mBurst End times: \033[0m\n" + str(array(time)[act_end]) + " s")

Burst Begin times: 
[ 1.56405484  4.87517094  8.26228971 11.80141381 14.7815183  17.40661035
 20.42271611 23.42282131 26.73593748] s
Burst End times: 
[ 2.44808584  5.60219644  8.86331079 12.5254392  15.49954348 18.42364601
 21.56775626 24.7688685  27.75497321] s

This procedure can be automatically done by detect_emg_activations function in detect module of biosignalsnotebooks package

In [24]: 
activation_data = bsnb.detect_emg_activations(signal, sr, smooth_level=20, threshold_level=10, time_units=True, volts=False, resolution=None, device=device, plot_result=True)

As described on the intro, electromyographic (EMG) signals are generated through voluntary actions of the subject, in contrast with electrocardiographic signals.

So, due to the voluntary nature, EMG signal is not being formed uninterruptedly and between muscular activations there are inactivation periods, consisting mostly in noise, which we want to avoid during our analysis.

With the steps described on the current Jupyter Notebook , user will be in possession of an important tool to start his EMG analysis.

We hope that you have enjoyed this guide. biosignalsnotebooks is an environment in continuous expansion, so don"t stop your journey and learn more with the remaining Notebooks !

☌ Project Presentation
☌ GitHub Repository
☌ How to install biosignalsnotebooks Python package ?
☌ Signal Library 		☌ Notebook Categories
☌ Notebooks by Difficulty
☌ Notebooks by Signal Type
☌ Notebooks by Tag
In [25]: 
from biosignalsnotebooks.__notebook_support__ import css_style_apply
css_style_apply()

.................... CSS Style Applied to Jupyter Notebook .........................
Out[25]: